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ABSTRACT 

Most field stars will have encountered the highest stellar density and hence the largest number of 
interactions in their birth environment. Yet the stellar dynamics during this crucial phase are poorly 
understood. Here we analyze the radial velocities measured for 152 out of 380 observed stars in the 
2-6 Myr old star cluster IC 348 as part of the SDSS-HI APOGEE. The radial velocity distribution of 
these stars is fitted with one or two Gaussians, convolved with the measurement uncertainties including 
binary orbital motions. Including a second Gaussian improves the fit; the high-velocity outliers that 
are best fit by this second component may either (1) be contaminants from the nearby Perseus OB2 
association, (2) be a halo of ejected or dispersing stars from IG 348, or (3) reflect that IG 348 has not 
relaxed to a Gaussian velocity distribution. We measure a velocity dispersion for IC 348 of 0.72 ±0.07 
km s“* (or 0.64 ± 0.08 km s“* if two Gaussians are fitted), which implies a supervirial state, unless 
the gas contributes more to the gravitational potential than expected. No evidence is found for a 
dependence of this velocity dispersion on distance from the cluster center or stellar mass. We also find 
that stars with lower extinction (in the front of the cloud) tend to be redshifted compared with stars 
with somewhat higher extinction (towards the back of the cloud). This data suggests that the stars 
in IG 348 are converging along the line of sight. We show that this correlation between radial velocity 
and extinction is unlikely to be spuriously caused by the small cluster rotation of 0.024 ± 0.013 km 
s“* arcmin“* or by correlations between the radial velocities of neighboring stars. This signature, if 
confirmed, will be the first detection of line-of-sight convergence in a star cluster. Possible scenarios 
for reconciling this convergence with IG 348’s observed supervirial state include: a) the cluster is 
fluctuating around a new virial equilibrium after a recent disruption due to gas expulsion or a merger 
event, or b) the population we identify as IG 348 results from the chance alignment of two sub-clusters 
converging along the line of sight. Additional measurements of tangential and radial velocities in IG 
348 will be important for clarifying the dynamics of this region, and informing models of the formation 
and evolution of star clusters. The radial velocities analyzed in this paper have been made available 
online. 

Subject headings: techniques: radial velocities open clusters and associations: individual IG 348 stars: 
pre-main sequence 


1. INTRODUCTION 

The early (i.e., few Myr) dynamical evolution of star 
clusters is still poorly constrained. Prestellar and pro- 
tostellar cores have a small velocity dispersion of 300 - 
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500 m ( Andre et al.|2007 Kirk et al.|2007 ), typically 
smaller than the dispersion measured from the linewidths 
of low-density gas tracers. Kirk et al. (2007) suggested 
that this dispersion is smaller than the dispersion needed 
to prevent collapse. This subvirial initial state should 
cause stars formed from the prestellar cores to fall into 
the potential well of the molecular cloud after they decou¬ 
ple from any support of the surrounding gas (e.g. from 
the magnetic fields). Indeed young stars typically a ppear 
to have a larger velocity dispersion of 1-2 km s“* dJoer- 
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gens 


20061 Furesz et al. 


however tnese stellar ve. 


2006 2008 Tobin et al. 2005) 


ocities have not been measured 


in the same regions as the velocities of the pre-stellar 
cores, precluding a direc t comparison. In a companion 
paper (Foster et al.|2015 ) analyze the radial velocity dis- 
tribution of JNGG 1333 to show that the stellar velocity 
dispersion is larger than the velocity dispersion of the 
pre-stellar cores in the same region and hence provide 
evidence for an initial global collapse (although there are 
alternative scenarios to explain the increase in the veloc¬ 
ity dispersion, such as the increase of the velocity disper¬ 
sion from the dynamic ejections from unstable multiple 
systems). 
























2 


If such an initial collapse takes place it should only 
last for about a free-fall time. Afterwards the dynamical 
evolution of the cluster will be dominated by mass loss, 
both stellar mass loss from for example stellar winds and 
gas mass loss from the dissipation of the natal molec¬ 
ular cloud. The latter is referred to as gas expulsion 
and is traditionally thought to play a crucial role in the 
evolution of a young embedded cluster, because in most 
star-forming regions up to ~ 10% of the molecular gas 
gets co nverted into stars (i.e., overall star-format ion effi¬ 
ciency; Evans et al. 2009). Simulations and analytical ap- 
proximations have shown that if the stars have the same 
spatial distribution as the gas and are in virial equilib¬ 
rium the star-formation efficiency will ha ve to be greater 


than 20-40% to survive gas expulsion (Tutukov 


1978 


Hills||1980[ Elmegreen||198^ |Lada et ah] 1984[ [Goodwin] 


& llastian|2006p. So local clusters are unlikely to survive 
gas expulsion, unless either (1) their star-formation effi¬ 
ciency is much higher than average, (2) they are still sub- 
virial at the time of gas expulsion (Goodwin 2009| , (3) 
or the stellar group has contracted into the center of the 
molecular cloud, cre ating locally a higher effective star- 


forma tion efficiency (Smith et al. 2011 Kruijssen et al. 

ml. 


To observationally constrain the dynamical evolution 
during this crucial epoch, we need to observe stellar ve¬ 
locities in clusters before and shortly after this gas ex¬ 
pulsion. Here we analyze stellar radial velocities derived 
as part of the INfrared Spectra of Young Nebulous Clus¬ 
ters (IN-SYNC) ancillary program of the Apache Point 
Observatory Galactic Evolution Experiment (APOGEE; 
Zasowski et al.|[2013 and Majewski et al. in p rep.) from 


the third S l oan Digital Sky Sur vey (SDSS-III; Eisenstein 
et al.|20lT ). Foster et al. (2015) analyzed the velocity dis¬ 


tribution in MGG 1333 and Ua Rio et al. (in prep) will 
look at the velocity distribution in Orion A. Here we fo¬ 
cus on IC 348, a young star-formi ng region in th e Perseus 
molecular cloud (see review from Herbst 2008). Age es- 
timates for the stars in IC 348 va ry from 2-3 Myr (e.g., 
Luhman et al. [2003) to ^ 6 Myr (Bell et al. 2013). For 
the velocity dispersion and half-rnass radius presented 
in this work we find that the stars cross the half-mass 
radius in 0.7 Myr, so for both age estimates IC 348 is 
many crossing times old. This implies th at any initial 


spati al substructure should have dissipated (Parker et al. 


stellar distribution was found 

in IC 348 (iMuench et al. 

2007 

Kumar & Schmeja|2007 

Schmeja et al.|2o011|). Any 


above should also have ended on the free-fall timescale 
of about 1 Myr and the subsequent evolution should be 
dominated by mass loss from winds and perhaps super¬ 
novae (altho ugh no evidence of supernova bubbles have 


been found; Ridge et al. 2006b) and the dissipation of 
the molecular cloud out ot which IC 348 formed. This 
implies that IC 348 should either be in virial equilibrium 
or should be expanding due to recent mass loss. In this 
paper we will show that the velocity dispersion of IC 
348 does indeed suggest that the cluster is supervirial. 
However, we will also show that the stars in IC 348 are 
actually converging along the line of sight, despite its ob¬ 
served virial state and the theoretical arguments above. 

In Section [2] we will discuss the selection of the 380 
observed stars, as well as the subset of 152 stars for 


which we analyze the radial velocity distribution. Here 
we also summarize the observations and the derivation of 
the stellar parameters including the radial velocity from 
high-resolution H-band APOGE E spectra. The spectra l 
analysis is described in detail in Cottaar et al. (2014). 
In Section P we discuss how we analyze the observed ve- 
locity distnbution in order to derive the mean velocity 
and velocity dispersion, corrected for measurement un¬ 
certainties and the velocity offsets due to binary orbital 
motions. The resulting mean velocity and velocity dis¬ 
persion and their dependence on other stellar parameters 
such as stellar position, mass, and interstellar extinction 
are discussed in Section (for the velocity dispersi on) 
and in Section (for the mean velocity). In Section |5.1| 
we will present the main result of the paper, namely that 
stars with larger extinctions are on average blueshifted 
compared with stars with lower extinction. In Section [^ 
we will look at possible scenarios to explain this correla¬ 
tion and conclude that the stars in IC 348 are probably 
converging along the line of sight. There we also discuss 
the virial state of IC 348 and interpret a velocity gra¬ 
dient observed across the plane of the sky. Finally our 
conclusions are presented in Section 

2. OBSERVATIONS 
2.1. Target selection 

Here we briefly discuss how the targets for the 
APOGEE/IN-SYNC spectroscopic survey of IC 348 were 
selected. We first describe the creation of the target cata¬ 
logue and analyze its completeness. Then we discuss the 
completeness of the subset of the stars in the catalogue 
for which spectra were actually taken. 

Likely IC 348/Perseus members were selected from a 
catalog of confirmed and candidate IC 348 members com¬ 
piled by August Muench (private communication), sup¬ 
plemented with mid-infrared excess sources id entified by 
the cores-to-discs (c2d) Spitzer survey team (Jprgensen 


et al.j 2006 Rebull et al. 2007). The Muench catalog 


includes 449 confirmed and candidate IC 348 members: 
the majority of these sources were drawn from previously 
published catalogs of photometrically selected, and of¬ 
ten spectroscopically confirmed, members (e.g 


Luhman 


et alj|1998[ Luhman||1999| Luhman et al. 2005 |21){M 
Muench et al.|2003 2007). Including candidates selectea 


on the bas is ot deep N-ray observations of the cluster 


cente r (i.e. Preibisch & Zinnecker 2001 
2003), this catalog is likely highly comp^ 


Preibisch et al. 
ete m the clus¬ 
ter ce nter (more than 80% complete for 77 < 16 accord¬ 
ing to Muench et al. 2007), as verified by our indepen¬ 
dent analysis ot source counts in the region (see Figure 
2, below). In the more poorly studied outskirts of IC 
348, the Muench catalog does include additional candi¬ 
dates selected via an R vs. R-J color-magnitude cut per¬ 
formed utilizing photometry from the USNO NOMAD 
catalog; these candidates are less secure, particularly at 
early types, where member selection via the R vs. R-J 
CMD becomes less efficient due to the steep cluster se¬ 
quence, and further from the cluster core, where source 
contamination from background stars and members of 
the Perseus OB2 association is expected to become sig¬ 
nificant. 

We provide evidence for the completeness of the target 
catalogue in the cluster core by studying the distribu- 
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tion of the 2MASS sources tha t are not in our targe t 
catalogue in a similar manner as Cambresy et al. (2006). 
If the target catalogue is complete, all stars not in the 
catalogue should be background and hence should have 
the same brightness distribution and spatial density as 
background stars in an off-cluster field (after correcting 
for extinction). The upper panel in Figureshows the 
extinction-corrected H-band density distribution of these 
untargeted 2MASS sources in the cluster core (blue) 
and the 2MASS sources in an off-cluster field covering 
20’ to 60’ from the cluster center (cyan). In the mag¬ 
nitude range covered by our spectroscopic observations 
{H < 13.5), the distributions are consistent, suggesting 
that the targeting catalogue is (mostly) complete. This 
is confirmed in the spatial distribution of non-targeted 
2MASS sources with H < 13.5 (lower panel of Figure!^, 
which shows no remaining overdensity of sources at me 
location of IC 348. 




RA (deg) 


Fig. 1.— The upper panel shows the stellar density distribu¬ 
tion per bin of H-band, corrected for the background extinction 
as estimated from the submillimeter continuum observed by Her- 
schel. The various histograms show the stellar density of the stars 
in the target catalogue (green), the stars actually observed (red), 
the 2MASS sources classified as background within the cluster cen¬ 
ter (within 6.6 arcminutes; blue), and 2MASS sources outside of 
the cluster (further than 20’; cyan). The bottom panels shows the 
density distribution of all 2MASS sources not classified as poten¬ 
tial members with H< 13.5. The two circles are at 6.6 (dotted) 
and 20 (dashed) arcminutes from the cluster center (same as in 
Figure!^. The lack of an overdensity of 2MASS sources classified 
as background at the position of IC 348 provides evidence for the 
completeness of the adopted targeting catalog. 


Although the target catalogue appears to be reason¬ 
ably complete (at least in the cluster core), the spectro¬ 
scopically observed stars are not complete. Most targets 


in the outskirts were observed, however only about half of 
the stars in the t^et catalogue in the cluster core were 
observed (Figure [^, because APOGEE is unable to si¬ 
multaneously target stars within ~ 71.5” of each other 
due to a collision through fiber assignment. This meant 
that even though 9 separate plates were drilled to cover 
IC 348, many stars in the dense cluster core could still 
not be targeted. 
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Fig. 2.— Map of the stars in the target catalogue in green. Stars 
circled in cyan have been observed. Additional background stars 
observed to fill up the APOGEE fibers have not been marked. The 
lower panel shows a zoom-in of the center of the upper panel, il¬ 
lustrating the lower coverage in the cluster center. For comparison 
with FigurelUwe added two circles at 6.6 (dotted) and 20 (dashed) 
arcminutes from the cluster center. T he bac kground shows a Her- 
schel column density map (see Section |4.3.2[ |. 

The priority of the targets was based on the H-band 
magnitude. The highest priority was given to targets 
with 7 < H < 12.5. For these high priority target we 
are mostly complete (^ 90%), even in the cluster center 
(upper panel in Figure [^. For fainter targets the prior¬ 
ity was assigned based on the H-band magnitude after 
correction of the background extinction, where the ex¬ 
tinction was estimated from the J-H vs. H-Ks 2MASS 
color-color diagram. Intrinsically brighter stars were tar¬ 
geted first. 

Most (90 %) stars were observed for multiple epochs. 
Half of these were observed within a single year (with a 
baseline up to a few months), while the other half were 
observed over two years with a baseline up to 500 days. 
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The spectra have a median signal to noise of 70. 


2.2. Spectral analysis 

Cottaar et al. (2014|) describes the analysis of the 
high-resolution near-intrared spectr a obtained by the 


APOGEE multi-object spectrograph (jWilson et al.(2012 1 
from stars in IC 348, NGC 1333, NGC 2264, and Orion 
A as part of the IN-SYNG an cillary program. In sum¬ 
mary BT-Settl model spectra (lA llard et al.||2011|) w ere 
fitted to every reduced spectruinl INidever et al.|2U15 1 in 
a minimum chi-squared sense after masking any bad pix- 
els or strong telluric emission lines. From our best-fit we 
extract the following parameters: the effective tempera¬ 
ture, surface gravity, H-band veiling, rotational velocity, 
and radial velocity. During these fits the observed con¬ 
tinuum was matched to the continua of the model spectra 
using a polynomial. Accurate estimates of the interstellar 
reddening were obtained by computing the E(J-H) of the 
stars in IC 348 with respect to an empirical color locus 
in the Pleiades, which we mapped to our effective tem¬ 
perature scale by measuring effective temperatures from 
APOGEE spectra observed for sta rs in the Pleiades us- 
ing the same pipeline as for IC 348. Cottaar et al. (20141 
showed for a subset of the stars in IC 348 tor which opti¬ 
cal photometry was available that the E(J-H) measured 
in the infrared accurately predicted the reddening ob¬ 
served in the optical, which put an upper limit on the 
reddening uncertainty of only about 0.06 in E(J-H). 

Throughout th is paper we wil l focus on the observed 
radial velocities. Cottaar et al. (2014) found a system¬ 
atic redshift for the coolest stars in IC 348 (as well as 
NGC 1333) of a few km s“^, which we argued was likely 
caused by inaccuracies in the molecular (probably wa¬ 
ter) line lists affecting stars cooler than 3300 K. An 
empirical spline was fitted to this systematic offset and 
subtracted to get a consistent zero-point across all stel¬ 
lar masses. This corrected radial velocity will be used 
throughout this paper. The radial velocities were not 
corrected for variations of the velocity zero-point over 
time, which were found to be small (up to about 100 m 
s—1) in the APOGEE survey (Nidever et al.||2015). The 
current APOGEE pipeline estimates typical radial veloc- 
ity uncertainties for our sample of around 100 m s“^. In 
our analysis, we scale those uncertainties to match the 
sample’s measured epoch-to-epoch RV variability, which 
is larger by a factor of ^3 than the APO GEE pipeline’s 
estim ates (see discussion of this effect by (Cottaar et al. 
2014) in their Sec. 3.2.1). As a result, the kinematic 
results that follow should be robust against all random 
uncertainties. The radial velocities (Table have been 
made available online with this publicatio n. The other 
parameter have been previously released by|Cottaar et al. 
|M4|. 


2.3. Sample for velocity analysis 

We make several cuts to the total population of ob¬ 
served cluster members to define a population of IC 348 
members useful for a dynamical analysis of the cluster. 
The main goal of these cuts is to only analyze stars for 
which the radial velocity has been accurately determined 
and for which this radial velocity is likely to be affected 
by the gravitational field of IC 348. First we make two 
cuts to the individually observed spectra: 


1. From the 2323 observed spectra in IC 348 we cut 
591 spectra with S/N less than 20. These spectra 
are so noisy that they provide very little informa¬ 
tion on the actual radial velocity of the star. Fur¬ 
thermore the convergence in the spectral fitting is 
often very poor for these low S/N spectra, which 
leads to offsets in the stellar parameters larger than 
expected from the already large noise level. 

2. From the remaining 1732 spectra we cut an ad¬ 
ditional 17 spectra with best-fit effective tempera¬ 
ture less than 2400 K. Our survey does not go deep 
enough to actually observe such cool stars, which 
suggests that these rare low best-fit effective tem¬ 
peratures represent fits to the noise in the spectra 
and not fits to the spectral features. Indeed all of 
the stars from which these spectra were taken have 
other epochs with more reasonable stellar parame¬ 
ters including larger effective temperatures. 

For 67 out of 380 stars observed in IC 348, no spectra 
survive the cuts described above. So after these initial 
cuts we are left with 313 stars. This is further reduced 
to a total used sample size of 152 by cutting any stars 
that fall within one or more of the following categories: 

1. 12 stars with vsini < 5 km s“^ are cut, because 
they are very unlikely to be actual cluster mem¬ 
bers. These stars are plotted in purple in Figure 

They have a velocity distribution with no clear 
peak around the velocity of IC 348 (top panel) and 
are typically hot stars (KS-test p-value = 10“® that 
they have the same effective temperature distribu¬ 
tion as the other stars), which are relatively faint 
(bottom panel). These stars are likely contami¬ 
nants, which ended up being observed based on 
their position in the R vs. R-J CMD. 

2. 11 stars with vsini > 150 km s“^ are cut. These 
are likely to be cluster members, however most of 
these stars appear to be very early type with 8 out 
of 11 stars having a best-fit effective temperature 
of 7000 K (black points in Figure!^, which is the 
upper bound in the adopted grid ofmodel spectra. 
The predominance of the hydrogen lines for these 
early-type stars might have caused an overestimate 
of their large rotational velocities, as well as a lower 
reliability of the radial velocities of these stars as 
shown by the increased spread in RV for these stars 
in the upper panel of Figure This increased un¬ 
certainty in the radial velocities is the main reason 
to exclude these stars. 

3. We cut 63 stars, which are more than 20’ (> 3 
half-mass radii) from the cluster center, because 
they are far enough that they are unlikely to be 
gravitationally affected by the cluster and hence 
do not provide us information about the current 
dynamical state of the cluster. Many of these stars 
show clear signs of youth and they might represent 
a population from distributed star-format ion in the 
Perseus molecular cloud, they might be related to 
the Perseus OB2 association in which IC 348 is 
embedded, or they might be ejected by IC 348. We 
will briefly discuss these outlying stars in Section 

Sjl 
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TABLE 1 

Sample of the radial velocities of stars in IC 348. The fourth colu mn is the radial velocities corrected for the 

SYSTEMATIC REDSHIFT FOUND FOR THE COOLEST STARS ||COTTAAR ET AL.|2014||. ThE FULL TABLE HAS BEEN MADE AVAILABLE ON THE 

JOURNAL WEBSITE AND HTTP : //WWW . ASTRO . UFL . EDU/lNSYNC/ 


2MASS 

Julian date 

radial velocity (km/s) 

corrected Urad (km/s) 

uncertainty Vrad (km/s) 

2M03233718-I-3056336 

2456564.00 

-98.50 

-98.50 

0.27 

2M03233718-I-3056336 

2456607.75 

-98.33 

-98.33 

0.26 

2M03233787-I-3131094 

2456561.75 

-2.03 

-2.03 

0.31 

2M03233787-I-3131094 

2456236.75 

-1.95 

-1.95 

0.26 

2M03233787-I-3131094 

2456315.75 

-1.99 

-1.99 

0.24 

2M03234277-I-3053142 

2456564.00 

-20.57 

-20.57 

0.30 

2M03234277-I-3053142 

2456607.75 

-20.73 

-20.73 

0.12 

2M03234517-I-3109561 

2456674.50 

37.55 

37.55 

0.16 

2M03234831-I-3121526 

2456674.50 

25.74 

25.74 

0.20 

2M03235905-I-3101512 

2456671.50 

48.43 

48.43 

0.32 

2M03235905-I-3101512 

2456674.50 

48.11 

48.11 

0.28 

2M03235905-I-3101512 

2456564.00 

47.87 

47.87 

0.39 

2M03235905-I-3101512 

2456607.75 

47.91 

47.91 

0.36 

2M03240225-I-3103502 

2456674.50 

-36.80 

-36.72 

0.20 


4. We cut 89 stars with uncertainties on the radial ve¬ 
locity averaged over all epochs larger than 500 m 
s“^. These stars have velocity uncertainties com¬ 
parable to the velocity dispersion itself and hence 
their velocity distribution is not dominated by the 
intrinsic velocity distribution of IC 348, but by the 
measurement uncertainties. These measurement 
uncertainties are the mselves uncertain an d appear 
to be non-Gaussian dCottaar et al.|[2014[ ). As we 
will show in Section l4.2l the effect ot this cut on the 
fit to the velocity distribution is very small, because 
the measurement uncertainties are explicitly taken 
into account in the fitting procedure. 


5. Out of the remaining 172 stars after the cuts de¬ 
scribed above, 20 more are cut because they are 
RV-variable. To define a star as RV-variable we 
first compute the for a model where the RV is 
the same over all epochs: 


2 _ (RVi — 

r * 


( 1 ) 


where we sum the square of the radial velocity 
offset from the weighted mean (RVi — fi) normal¬ 
ized by the measurement uncertainty Ui over all 
epochs i. If the probability of drawing the com¬ 
puted or ^ larger due to random noise is 
smaller than 10“"^ we conclude the star is RV- 
variable. This probability is computed from a y^- 
distribution with the degrees of freedom set to the 
number of epochs minus one. 

None of these cuts are designed to match the mass or 
spatial distribution of the adopted sample to the actual 
mass or spatial distribution of all IC 348 cluster mem¬ 
bers. This means that if the stellar velocity dispersion 
depends on stellar mass or location, the velocity disper¬ 
sion of the sample analyzed here will not be representa¬ 
tive of the dynamical state of all stars in IC 348. We will 
show t hat these possible biases have a limit ed e ffect in 
section [T4] for the spatial bias and in section [T5| for the 
bias in stellar mass. 


3. METHOD TO FIT THE VELOCITY DISTRIBUTION 

The observed radial velocity distribution of our ana¬ 
lyzed sample of 152 stars (with a total of 758 observed 


i 


remaining stars 
4(7 RV outliers 


vsini <5 km S 
■ysini >150 km S“ 


-50 0 50 100 150 

mean radial velocity (km s ) 



Fig. 3.— Two plots of the stellar parameter distribution for slow 
rotators (purple; usinf < 5 km s“^) and rapid rotators (black; 
usin 2 > 150 km s~^). The remaining population of intermediate 
rotators is split into those stars with radial velocities within 3 km 
s~^ (i.e., ~ 4 times the velocity dispersion of IC 348) of the mean 
velocity (blue) and those with a radial velocity offset from the mean 
larger than 3 km s“^ (red). The top panel show the weighted 
mean radial velocity over all epochs versus the rotational velocity 
usinL The bottom panel shows the effective temperature versus 
the extinction-corrected J-band magnitude. 


spectra) in IC 348 is not only influenced by the stellar 
motions through the cluster’s potential, but also by or¬ 
bital motions of binary members and potential contami¬ 
nation of background or foreground stars, especially from 
young stars from the surrounding Perseus OB associa¬ 
tion. All these effects create a complex non-Gaussian ve- 
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locity distribution. We model this distribution to retrieve 
the underlying dynamical properties of IC 348 using the 
maxim um- likelihood procedure outlined in Cottaar et al. 
(2012) and Cottaar & Henault-Brunet (2014|). 

in summary we start by assuming that iC 348 has in¬ 
trinsically a Gaussian velocity distribution: 


ingly single stars 


-P(^'intr) = 






( 2 ) 


where P(uintr) is the probability distribution of the in¬ 
trinsic velocity of the star (uintr) and /ic and CTc are re¬ 
spectively the mean velocity and velocity dispersion of 
the cluster. To actually compare this with the observed 
velocity distribution we have to convolve this distribution 
of intrinsic velocities with the probability distribution of 
the measurement uncertainties (P(uunc)) and the binary 
orbital motions {P{vhin))- 

P { Vohs ) =(1 - /bin)-f’(^intr) * C’(l'unc) 

+ fhinPi'Vintr) * -P(Wunc) * P{vhin), (3) 

where is the binary fraction after removing RV- 
variable stars, which we will define in equation]^ This 
convolution will have to be computed separately Iot every 
star, because the measurement uncertainties and the ef¬ 
fects of binary orbital motions differ per star (see below). 
Because we only analyze the stars for which the radial 
velocity has been measured with a precision smaller than 
500 m s“^ the width of the measurement uncertainties 
{P{vunc)) is for all stars smaller than that of the intrinsic 
velocity distribution (P(uintr)), however for completeness 
we still do include it as a Gaussian distribution, which is 
only a rough approximation of the true probab ility dis¬ 
tribution d ue to the measurement uncertainties (Cottaar 

eTalllMIt . 


Computing the velocity offsets due to binary orbital 
motions (P(ubin)) is more complicated, because the mul¬ 
tiple epochs over which APOGEE spectra were taken 
allows the binaries with short periods (up to several 
times the observational baseline) to be detected through 
their radial velocity variations. However, many binaries 
with orbital velocity amplitudes comparable to or larger 
than the velocity dispersion of IC 348 remain undetected, 
broadening the observed radial velocity distribution. We 
already excluded from the analysis all stars (20 out of 
172) with variable radial velocities (i.e. inconsistent with 
being single at the p < 10“^ level according to a y^-test). 
For the remaining stars we compute the distribution of 
radial velocity offsets due to binary orbital motions by 
(1) drawing in a Monte-Carlo like fashion a large num¬ 
ber of binary orbits, (2) creating for every star artifi¬ 
cial RV observations given the observational epochs and 
measurement uncertainties for that star, (3) discarding 
those binary orbits that would have been detected in the 
artificial RV observations by the y^-test, and (4) com¬ 
puting the distribution of radial velocity offsets between 
the obser ved and systematic velocities for t he remaining 
binaries (Cottaar & H&ault-Brunet 2014). This radial 
velocity offset distribution is used as i-’(ubin) in the con¬ 
volution in equation]^ This same population of binaries 
is also used to estimate the binary fraction of the seem- 


fhin 


1 - /det 
1 /*bin/*det 


( 4 ) 


where /bin is the binary fraction for the primary mass of 
the star prior to any multi-epoch observations and /det is 
the fraction of potential binary orbits, which could have 
been detected in the multi-epoch observations (i.e., the 
fraction of simulated binary orbits discarded in step 3). 

The first step of this procedure requires the drawing 
of many random binary orbits and hence we need to as¬ 
sume a period, mass ratio, and eccentricity distribution 
for the observed stars. We are mainly interested in the 
number of binaries with periods between tens and thou¬ 
sands of years, whose binary orbital motions cause ve¬ 
locity offsets on the order of the velocity dispersion of 
IC 348. The binary distribution in this period range has 
only been extensively cha racterized in the solar ne igh- 
borhood (see review from Duchene & Kraus 2013|), so 
we will use these results to inform our assunied orbital 
parameter distributions. First we assign to every star 
an age and a mass based on the Dartmouth tracks on 
their effective temperature and extinction-corrected ab¬ 
solute J- band magnitude (as suming a distance modulus 
of 6.98; Ripepi et al. 2014). The binary fraction and 
the mean and width of the semi-major axis distribution 
all decrease tow ards lower masses (e.g. Lada 2006 Jan- 
son et al. 2012). For brown dwarves w e set the binary 
fraction to 22% (Burgasser et al. 2007). Above the hy¬ 
drogen burning limit we let the binary fraction increase 
linearly, so that the binary f raction at one solar mass 
is 44% ( Raghavan et al.||2010 ). We implement the semi¬ 
major axis distribution as three log-normals, with a mean 
of log On = 0.86 and a width of criogo = 0.24 below 
0.2 Mq (iBurgasser et al.||2007|), with log on = 1-2 and 
ginir a = 0.8 between 0.2 Mq arid 0.8 Mq (Janson et al. 


2012), a nd with logon = 1-64 and oioga = 1.52 above 
0.8 Mq (iRaghava n et al.||2010|). We adopt a flat mass 
ratio distribution (Reg giani fc Meyer 2013 ) an d flat ec¬ 
centricity distribution ( Duchene &: Kraus||2013|) over the 
whole mass range. These assumptions about the binary 
properties based on the field stars are expected to only 
app roximate the true bin ary properties in IC 348, how¬ 
ever Cottaar et al. (2012) showed that the best-fit mean 
velocity and velocity dispersion are not very sensitive to 
the adopted binary properties. 

The efficiency with which binaries could be detected 
(as estimated by /det) varies widely across the sample. 
Out of the 152 stars whose velocities are analyze d in this 
paper (i.e. after the cuts described in section 2.3) 20 
have only been observed for a single epoch, 33 tor two 
epochs, while the remaining 99 stars have been observed 
for three up to thirteen epochs. For the stars with mul¬ 
tiple epochs, the total observational baseline varies with 
all observations for 75 stars taken within a single season, 
while for the remaining 57 stars we have data taken over 
two years. Finally the RV precision of the observations 
also varied between 60 and 500 m s“^ with a median of 
200 m s—1. The wide variety of binary orbital parame¬ 
ter distributions assumed for different mass ranges and 
the variety of observational constraints lead to a binary 
detection rate (/det) that ranges from 0 for stars with a 
single epoch and a few percent for stars below 0.2 Mq up 
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to ^ 20% for solar-type stars with a baseline of at least a 
year. However, even in this best-case scenario (i.e. solar 
mass stars with a baseline of ^ 1 year) we can not detect 
binaries with orbital velocities comparable to the veloc¬ 
ity dispersion which would have periods of about 10^ 
years. As a result, despite this careful analysis of these 
multi-epoch observations, we still expect the majority of 
relevant binaries to go undetected, so that changes in our 
binary detection efficiency will have only a limited effect 
except in the wings of the measure d velocity distribution 
(Cottaar & Henault-Brunet||2014|). 

i'his procedure allows us lo compute the likelihood of 
observing a given velocity (H(fobs)) from the convolution 
in equation 1^ given a set of variables parameterizing the 
intrinsic velocity distribution (P(nintr))- We consider the 
best-fit parameters to be those where the total likelihood 
of reproducing all observed velocities (uobs) is maximized 
(i.e., maximize £ = 11 ^obs(i'obs))- The uncertainties 
on these parameters will be computed by Markov chain 
Monte Carlo (MCMC) simulations throughout this pa¬ 
per. To compare the relative accuracy of various mod¬ 
els of the mean velocity and velocity dispersion profile 
we will use t he Bayesian Information Criterion (BIC; 
Schwarz|1978 1. This criterion states that the best model 


is that, which minimizes 


BIC = —21n£ -I- khi{n), 


(5) 


where £ is the best-fit likelihood of reproducing all ob¬ 
served velocities, k is the number of free parameters, and 
n is the number of fitted datapoints (i.e., 152 radial ve¬ 
locities). The first term decreases if the fit is improved, 
while the second term acts as a penalty on the number of 
parameters, which prevents overfitting. The BIC values 
for the models discussed in this paper are shown in Table 

m 

Through the paper we vary on this general scheme in 
two important ways: we either take the mean velocity 
(/Tc) and the velocity dispersion (tJc) of IC 348 as global 
variables that are constant across IC 348 or we let them 
vary as a function of stellar parameters (e.g., stellar mass 
or location). The rows in Table represent different 
models for the mean velocity and/or velocity dispersion, 
which we shall discuss throughout the rest of this paper. 
The second important variation is whether we consider 
the possibility that a single Gaussian distribution might 
be insufficient to represent the true intrinsic velocity dis¬ 
tribution. To explore this we also model the intrinsic ve¬ 
locity distribution as a sum of two Gaussian, one of which 
has the IG 348 mean velocity (^c) and velocity disper¬ 
sion (uc) and the other of which has an alternati ve mean 
velocity (/Z g) and velocity dispersion (cto) (e.g., Jeffries 
et al.|[Ml4). Such a model has three more free parame¬ 


ters, namely the alternative mean velocity and velocity 
dispersion and the fraction of analyzed stars, whose ve¬ 
locity distribution is described by this second Gaussian. 
The columns in Table show the result for both a single- 
Gaussian model and a double-Gaussian model. 

4. TRENDS IN VELOCITY DISPERSION 
4.1. Global velocity dispersion 

We start with our simplest model fit to the observed 
velocity distribution, namely an intrinsic Gaussian ve¬ 
locity distribution with a constant mean velocity and 


velocity dispersion throughout the cluster, which is con¬ 
volved with both the measurement uncertainties and the 
binary orbital motions as described in section 1^ This 
model only has two free parameters, namely the mean 
velocity and velocity dispersion of IC 348. The resulting 
best-fit intrinsic Gaussian velocity distribution has been 
plotted in green in Figure 15 The posterior probability 
distribution of the mean velocity and velocity dispersion 
resulting from the MCMC simulations have been shown 
in green in the other two panels. Marginalized over the 
other parameters the MCMC simulation gives a mean he¬ 
liocentric velocity of 15.37±0.07 km s~^pqand a velocity 
dispersion 0.72 ± 0.06 km s“^. 

The observed velocity distribution does not appear to 
be truly Gaussian, however, due to an excess of stars at 
high and low velocities relative to the cluster’s mean ve¬ 
locity, compared to the number expected from a Gaus¬ 
sian distribution even after taking into account binary 
orbital motions. Although visually the small number of 
’’outliers” has a minor effect (left panel in Figure 4), the 
likelihood of reproducing all observed velocities increases 
significantly when the fit allows for a profile with stronger 
high velocity wings. This is true for tests of profiles with 
intrinsically stronger wings (i.e., a Lorentzian profile), 
but the high velocity outliers in our dataset are best 
modeled with a separate Gaussian component, whose 
amplitude can be tuned to the number of high velocity 
objects. Adding this second Gaussian leads to a signifi¬ 
cant decrease of 125 in the BIC (upper row in Table 
compared with the single-Gaussian model, which is muai 
larger than 10, the limit usually adop ted as providing 
decisi ve evidence in favor of a model (Kass & Raftery 
19951. This second Gaussian is found to have a mean 
velocity of 14 ± 4 km s“^ and a velocity dispersion of 
5 ± 4 km s“^ with only 10 ± 5% of the total stars con¬ 
tributing to this second Gaussian. The large parameter 
uncertainties are due to the small contribution of this 
Gaussian and due to its large overlap with the veloc¬ 
ity distribution of the first Gaussian. Adding this sec¬ 
ond Gaussian reduces the velocity dispersion of the main 
Gaussian to about 0.64 ± 0.08 km s“^ with no signifi¬ 
cant change in the mean velocity (Figure]^. Although 
the double-Gaussian model does provide a significantly 
better fit than the single-Gaussian model (according to 
the BIG), we will in this paper continue to consider both 
models, which generally agree very well with each other 
on the dynamical properties of 1C 348. 

The physical interpretation of this second Gaussian is 
unclear. It could possibly represent contamination from 
the Perseus OB2 association. However, the mean veloci¬ 
ties of the two Gaussians are consistent with each other, 
which suggests that they are a single population. In that 
scenario it could be that the second Gaussian represents 
stars ejected from IC 348 or that the velocity distribu¬ 
tion of IC 348 is simply non-Gaussian at this young age. 
We consider the latter option plausible, because we are 
unaware of any predictions that a turbulent molecular 
cloud would produce stars with a Gaussian velocity dis¬ 
tribution or that a non-Gaussian initial stellar velocity 

The uncertainty in the mean velocity is actually dominated 
by a possible systematic offset with resnect to RV s tandards on 


tire order of about 0.5 km s ^ (jCottaar et al. |2014l, not by the 
statistical uncertainty quoted here. 












TABLE 2 

Bayesian information criterion (BIC; eq. for the fits to the velocity distribution discussed throughout this paper. 


Description of fit 

Mean velocity 

Velocity dispersion 

Global fit 

Polynomial radial dispersion profile 
Polynomial dispersion profile with mass 
Step function at ~ 0.44 Mq 

Step function at ~ 0.8 Mq 

Velocity gradient with extinction 
Velocity gradient on the sky 

constant 

constant 

constant 

constant 

constant 
equation 
equation 14 

constant 

equation 
equation 
equation 12 
equation 12 
constant 

constant 


BIC (# parameters) 
Single Gaussian 
648 (2) 

658 (4) 

657 (4) 

656 (4) 

659 (4) 

636 (5) 

658 (4) 


BIC (# parameters) 
Double Gaussian 
523 (5) 

532 (7) 

531 (7) 

529 (7) 

542 (7) 

509 ( 8 ) 

532 (7) 


— best-fit single Gaussian — best-fit double Gaussians I I mean stellar velocity 




III_I_ 

0.4 0.5 0.6 0.7 0.8 0.9 1.0 


velocity dispersion (km s M 


Fig. 4.— Left panel: Histogram of the average velocities of all stars included in the fit (blue) with the best-fit intrinsic velocity distribution 
from the single-Gaussian model (green) and the double-Gaussian model (red). The effect of binary orbital motions on the best-fit models 
has been included in the fits, but not in the plot, because this effect is different for every star depending on the precision of the velocity 
measurements and the baseline over which they were observed. Other panels: The probability distribution in the mean velocity and 
velocity dispersion for the single-Gaussian (green) and the main peak of the double-Gaussian (red) model. The vertical colored lines in 
the right panel mark the velocity dispersion expected for virial equilibrium (dashed colored lines) and a dynamically unbound state with 
a net positive energy (towards the right of the solid colored lines) under various assumptions for the contribution of the gas mass to the 
dynamical state of IC 348, namely yellow for no contribution from the gas, cyan for a small contribution of 80Mq, and gray for a large 
contribution of 21 OM 0 (see Section|4.3.3[|. 


distribution would relax to a Gaussian one within the 
few Myrs, as required given IC 348’s age. Recent expan¬ 
sion or collapse of IC 348 might also have altered the 
velocity distribution. Even though the velocity distribu¬ 
tion of IC 348 appears to be non-Gaussian, we still adopt 
the Gaussian model, because it gives us two easily inter¬ 
pretable values (i.e. the mean velocity and the average 
deviation from this velocity). For the single-Gaussian 
model the velocity dispersion can be interpreted as the 
average deviation from the mean velocity for all IC 348 
stars, while for the double-Gaussian model the velocity 
dispersion represents the average deviation for the subset 
of IC 348 stars with similar velocities, while ignoring the 
relative large contributions of a few stars with outlying 
radial velocities. 

In this section we will continue to discuss the velocity 
dispersion and what it means for the dynamics of IC 348. 
In section we will then focus on the mean velocity and 
its spatial dependence. 


4.2. Effects of the sample selection 

H ere w e look at the effect of the cuts described in Sec¬ 
tion 2.3 on the fit to the velocity distribution. In par¬ 
ticular we will focus on two cuts, that remove a large 
number of stars: the 89 stars cut because of their large 
RV measurement uncertainty (> 500 m s“^) and the 63 
stars cut because they lie more than 20' from the cluster 
center. 

Figure [5] illustrates the changes in the single-Gaussian 


fit if the stars with large RV measurements are included. 
Because the RV uncertainties are explicitly included in 
the fit, the inclusion of the high-uncertainty stars actu¬ 
ally has very little effect with shifts of ^ 20 m s“^ in 
the mean velocity and velocity dispersion and minimal 
reductions in the precision of these parameters. For the 
high-uncertainty stars the offset from the mean velocity 
is mostly set by the measurement uncertainty rather than 
the velocity dispersion, so they provide little to no infor¬ 
mation on the velocity dispersion. Similarly, their larger 
spread implies that they do not provide a lot of infor¬ 
mation about the mean velocity of the cluster. When 
a double-Gaussian fit is considered there is a significant 
difference. If the stars with large RV measurements are 
included, the best-fit second Gaussian has a contribution 
of 4 ± 2% with a mean velocity of 3 ± 20 km s“^ and a 
velocity dispersion of 32±43 km s“^ (data not shown), so 
this second Gaussian now only provides a fit to outlying 
(probably bad) radial velocities, rather than suggesting a 
non-Gaussian velocity distribution of IC 348 like above. 
Adding this second Gaussian does not change the mean 
velocity or velocity dispersion of the first Gaussian in a 
significant way. Because of these outliers and the minor 
contribution of these stars with large RV measurements 
on the fit, we retain our cut in the measurement uncer¬ 
tainties for all future fits. 

The results of fitting the RV distribution of the stars 
beyond 20' (after making the same quality cuts) have 
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— Gaussian fit to (Tjiv <0.5 km s ^ — Gaussian fit to all o-j^v I I (Tjiv<0-5kms ^ I I all o-j^v 



Fig. 5. — Left panel: Histogram of the velocities of all stars included in the fit with the cut in the RV measurement uncertainty at 500 m 
s“^ (blue) and without this cut (orange). The lines show the best-fit single-Gaussian intrinsic velocity distribution to the radial velocities 
with the measurement uncertainty cut (green) and without the cut (magenta). Middle panel: Posterior probability distribution of the mean 
velocity with the cut in the RV measurement uncertainty (green) and without the cut (magenta). Right panel: Same as the middle panel, 
but with the posterior probability distribution of the velocity dispersion rather than the mean velocity. 


been illustrated in Figure The RV distribution of 
these spatial outliers is poorly fitted by a single Gaus¬ 
sian (see green line in left panel) and hence requires a 
double-Gaussian fit. These Gaussians seem to have dif¬ 
ferent mean velocities (Figure [^, which suggests that 
they represent two different populations. One of these 
populations (41 ± 11% of stars; fi = 18.0 ± 1.3 km s“^; 
tr = 4.0 ± 1.1 km s“^) appears to be inconsistent with 
either the velocity distribution of IG 348 or the much 
wider field velocity distribution (see e.g. slow rotators in 
upper panel of Figure]^. The mean velocity is consistent 
with the mean velocity of 23.5 ± 3.9 km s ~^ previously 


measured f or the Perseus OB2 association (Steenbrugge 
et al.|2003|) , which covers the observed region on the sky. 
The other population (59 ± 11% of stars; /r = 15.5 ± 0.2 
km s“^; a = 0.4 ± 0.2 km s“^) is consistent with the RV 
distribution of IC 348, which suggests that either these 
stars are ejected from IC 348 or they formed outside of IC 
348 from the same molecular cloud. If this extended pop¬ 
ulation is also present in front and behind IC 348 along 
the line of sight, it might explain the non-Gaussianity 
observed in the radial velocity distribution in IC 348. 

4.3. Mass and virial state of IC 348 
4.3.1. Stellar mass 

We need an estimate of the stellar and gas mass of 
IC 348 to determine its virial state. First we derive the 
stellar mass in IC 348. We will only consider stars more 
massive than 0.25 Mq (i.e., hotter than ^ 3300 K), be¬ 
cause these are typically brighter than 11=13 for which 
our target catalogue of potential members see ms to be 
complete (see completeness analysis in Section 2.1 and 
Figure [^. 

First we derive the stellar mass of the stars with 
APOGEE spectra using the spectroscopic effective tem¬ 
perature and extinction-corrected absolute J-band mag¬ 
nitude from the Dartmouth isochrones. Here we include 
all stars targeted as potential members, except for the 
slow rotators (zisin* < 5 km s“^), which are probably 
non-members as we argue in Section |2.3[ Eor the stars 
with spectroscopic masses above 0.25 Mq the total mass 
adds up to about 145 ±32 Mq within 20”, where the un¬ 
certainty comes from an assumed system atic uncertainty 
of 200 K on the effective temperature (Cottaar et al. 


20141. 

We then extend this mass estimate to the unobserved 
stars in the target catalogue. We assume that these un¬ 
observed stars have the same properties as the observed 
stars at the same apparent H-band magnitude, because 
the priority in the targeting was based on this magnitude. 
For every unobserved star we select the 20 observed stars 
with the most similar apparent H-band magnitude and 
randomly assign one of their masses to the unobserved 
star. On average this leads to an additional mass of 24±5 
Mq for cluster members with a mass above the 0.25 Mq. 
Hence we estimate the total cluster mass for stars more 
massive than 0.25 Mq to be 169 ± 37 Mq. 

Finally we estimate how much of the total cluster mass 
we missed by assigning a minimum stellar mass of 0.25 
Mf7^. For this we use the initial mass function (IMF) from 
Maschberger (2013), whic h has be e n des igned to match 
the classica l IMl*”s from Kroupa (2001) and Chabrierf 
(2003 20051, but has the advantage of consisting of a sin¬ 
gle mathernatically convenient continuous function. By 
integrating this IM F with the default parameters from 
Maschberger (2013) we find that stars with masses be- 
low 0.25 Mq contribute 17% of the initial total mass. 
Correcting for this implies a final total stellar mass in IC 
348 of 204 ± 45 Mq. The distribution of this mass has 
been plotted in Figure]^ in red. 

4.3.2. Gas mass 

We estimate the enclosed gas mass in IC 348 as a func¬ 
tion of projected distance from the cluster center using 
three different methods. In the first, we use the publicly 
available Herschel data for Perseus from the Gould Belt 


Survey (Andre et al. 2010) which measures the thermal 
emission from dust at 100, 160, 250, 350 and 500 /rm. We 
convolve all these maps to the resolution of the coarsest 
map (500 p,m, 37”) and for each pixel fit for the tem¬ 
perature and column density using a greybody fit, where 
the opacity/emissivity at each wa velength is taken from 
the Ossenkopf & Henning (1994) model for dust with 
thin ice mantles coagulating for 10® years at a density 
of 10® cm“®. This differs from the st andard procedure 
of fitting submillimeter emission (e.g. Hildebrand 19831 
because rather than assuming that emission Iroiii dust 
is a black-body modified by some power-law with expo- 
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— best-fit single Gaussian — best-fit double Gaussians I I d < 20' I I d > 20' 



Fig. 6 .— Left panel: Histogram of all stars included in primary fit (i.e., within 20^ from cluster center; blue) and all stars with the same 
cuts, but beyond 20' from the cluster center (orange). The best-fit single-Gaussian (green) and double-Gaussian (red) model to the stars 
beyond 20' from the cluster center have been overplotted. Middle and right panel: Posterior probability distribution in the mean velocity 
and velocity dispersion of the single-Gaussian (green) and double-Gaussian model (red) for the stars beyond 20' from the cluster center. 
For the double-Gaussian model the red lines show the posterior probability distributions of both Gaussians (solid and dashed). 



Fig. 7.— The enclosed stellar and gas mass in IC 348 as a function 
of distance from the cluster center. For the gas mass estimates from 
the Herschel dust emission maps, m aps, and maps have 

been included. The core and halo radius of|Herbst|(|2008|l have been 
marked together with the half-mass radius adopted here of 6.3’. A 
dense gas clump called the southwestern ridge extends from about 
500 to 700” to the south of IC 348. 


depth of the (1-0) emission. These relationships are 
calibrated for different sub-regions in Perseus, and we use 
the parameters for IC 348 but note that the calibration 
in this region actually excludes the core of the cluster 
where the COMPLETE near-infrared extinction map is 
unreliable due to the large density of non-background 
stars. 

Figure shows these three different estimates for 
the enclosed mass, which show substantial disagreement 
among the methods, illustrating the fundamental uncer¬ 
tainty of this calculation. We adopt the Herschel dust 
mass primarily because the CO conversion relations are 
not well calibrated in the central region of the cluster 
where we most care about the enclosed gas mass. The 
presence of substantial prot ostellar outflows and sh ells in 
the vicinity of the IC 348 (Arce et al. 2010, 2011) may 
also be complicating the calculation of gas mass from CO 
emission. 


4.3.3. Virial state 


nent /3, it allows the dust emissivity to have the more 
complicated dependence on w avelength estimated by the 


Ossenkopf & Henning (1994) model. In addition, we do 
not assume that the emission is optically thin. We still 
assume a single temperature along the line of sight and 
that the dust optical properties are constant along the 


line of sigh t (and well described by the Ossenkopf & Hen-, 
T994| model), with all the concurrent unce rtainty 


nmg 


(e.g. IShetty et al. |2009|). Following ||Hildebrand| (|1983 ), 
we adopt a gas-to-dust mass ratio of 101 ). 

We compare two alternate estimates of the ga s mass 


in IC 348 us ing CO maps from the COMPLETE (Ridge 
et al.|20dfe) Survey. We use the relationships calibrated 


Pineda et al. (2008) to convert the CO emission to 


a total gas mass; these relationships are calibrated with 
reference to the COMPLETE near-infrared dust extinc¬ 
tion maps, and thus ultimately also represent a dust mass 
with an assumed gas-to-dust mass ratio of 100. We use 
both the non-linear relationship to convert integrated in¬ 
tensity of ^^CO (1-0) to gas mass (Eq. 19) and the lin¬ 
ear relationship between the column density of ^^CO ( 1 - 
0) and gas mass (Eq. 21), where the latter calculation 
assumes that the excitation temperature of ^^CO ( 1 - 0 ) 
and ^^CO ( 1 - 0 ) are equal in order to estimate the optical 


To determine whether IC 348 is expected to expand we 
compare the measured velocity dispersion of 0.72 ± 0.06 
km s“^ (or 0.64±0.08 km s“^ when fitting two Gaussians) 
to the velocity dispersion needed for virial equilibrium 
assuming spherical symmetry, which is given by: 


CTdyn 



( 6 ) 


where G is the gravitational constant, M is the total stel¬ 
lar and gas mass, rhm is the half-mass radius, and 77 is a 
structural parameter that depends on the d ensity distri¬ 
bution . To estimate rhm and rj we fit an |Elson et al.| 
(1987) profile to the stellar distribution. This profile 
has a flat surface density distribution in the cluster core, 
which transitions to a power-law surface-density profile 
with slope —7 at a distance a: 


S = So 



(7) 


where R is the projected distance from the cluster center 
and Eq is the central surface density. To minimize the 
extinction bias from the patchy molecular cloud on this 
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fit we only include the stars with extinction-corrected 
J-band magnitude brighter than 12.5. Because we only 
have an extinction-corrected J-band measured for stars 
with APOGEE spectra, this limits ourselves to the sub¬ 
sample of observed member stars. This should be a minor 
effect, because the vast majority of these stars have an 
apparent H-band magnitude brighter than 12.5, which is 
a magnitude range over which we are mostly complete 
(Figure [^. We only fit the profile out to 20' from the 
cluster center to minimize the effect of individual out¬ 
liers. The density profile is fitted to the data by maximiz¬ 
ing the likelihood (£) of reproducing all stellar distances 
from the cluster centers: 


marked as vertical lines in the right panel of Figure 15 
The observed velocity dispersion was 0.72 ± 0.07 Km 
s“^ for the single-Gaussian model and 0.64 ± 0.08 km 
s“^ for the double-Gaussian model. This is sufficiently 
high that the gas will have to significantly contribute to 
the gravitational potential for the cluster to be in virial 
equilibrium, with the observed velocity dispersion being 
more than 2 -sigma higher than the velocity dispersion 
estimates with no gas or only a little gas (see Table [^. 
This analysis suggests that IC 348 might be slightly su- 
pervirial, however the velocity dispersion does not appear 
to be so large that it can not regain virial equilibrium 
through expansion. 


C = AY[R,nA), ( 8 ) 

i 


where we multiply over all stars f at a distance Ri 
from the cluster center and the normalization of the 
probability distributions has been summarized by a sin¬ 
gle constant A. This fit results in a = 4.1’ and 
7 = 3.1, which corresponds to a half-mass radius of 
6.3’ or 0.47 pc at a distance to IG 348 of 273 ± 23 pc 
(Ripepi et al.||2014), where we assume that the distribu- 
tion of stars with extinction-corrected J-band magnitude 
brighter than 12.5 is representative of the distribution of 
all stars. This matches the observed half-mass radius 
(Figure reasonably well, which provides support to 
the adopted density prohle. iP ortegies Zwart et al .| (|2010| 
showed that for 7 > 2.5 the |Flson et al.j ( |1987| ) prohle 
consistently yielded rj ~ 10 , which we will adopt here. 

Entering these values in equation we hnd a velocity 
dispersion of 0.44 ± 0.06 km s“^ for virial equilibrium if 
we ignore the gas mass, while a cluster with a velocity 
dispersion higher (i.e. ~ 0.6 km s“^) would have a 
positive total energy, which could cause the cluster to un¬ 
bind. Adding the gas mass to this estimate of the virial 
state is more tricky, because the gas is clearly not spher¬ 
ically distributed around the center of the stellar cluster 
with a large fraction contained in the south-western ridge 
(see Figure [^, where its gravitational potential is less 
likely to contribute to holding IC 348 together, than if 
this gas was found in the center of IC 348. Projected onto 
the region of interest (within 20' from the cluster center), 
there is a total of 2 IOM 0 of gas according to the Herschel 
dust maps, however only 4 OM 0 of this gas is actually pro¬ 
jected within one half-mass radius of the cluster center. 
Based on the spread in extinction observed for the cluster 
members in IC 348 we conclude that the stars are embed¬ 
ded within at least this gas. So as a lower estimate for 
the effect of the gas on the dynamical state of IC 348 we 
assume that this 4OM0 has the same density profile as 
the stars for a total extra mass of 80Mq, which leads to 
a velocity dispersion of 0.51 ± 0.04 km s“^ with the clus¬ 
ter being unbound for a velocity dispersion greater than 
^ 0.72 km s“^. Finally as an extreme upper case we con¬ 
sider if all the gas of ~ 2 IOM 0 had the same distribution 
as the stars (even though this is clearly inconsistent with 
the observed surface density distribution). In that case 
the velocity dispersion in virial equilibrium would be sig¬ 
nificantly higher at 0.61 ± 0.04 km s“^ with the cluster 
being unbound at ~ 0.87 km s“^. All these estimates of 
the velocity dispersion in virial equilibrium for different 
gravitational influences of the interstellar gas have been 


4.4. Velocity dispersion profile 

In globular clusters the velocity dispersion is often 
found to drop with increasing dista nce from the clus¬ 
ter center (e.g. Watkins et al. 20131 as expected from 
the Jean’s equation. Such a trend would bias our veloc¬ 
ity dispersion, because our completeness is lower in the 
cluster center, where the survey efficiency was limited by 
the collision radius of the input hbers. Here we explore 
how the velocity dispersion varies with distance from the 
cluster center under the assumption that this trend can 
be approximated by a second-order polynomial: 


log(crc) =Po+ Pid + P2d‘^ 


(9) 


where d is the distance from the cluster center and Pq, 
Pi, and p 2 are the polynomial constants. We fit the poly¬ 
nomial to the logarithm of the dispersion to ensure that 
the dispersion remains positive. The mean velocity is still 
kept the same for all stars, which leads to a total of four 
parameters for the single-Gaussian model and 7 for the 
double-Gaussian model. Using this parameterization of 
the velocity dispersion and mean velocity we then maxi¬ 
mize the likelihood of reproducing the observed velocity 
distribution (equation]^ and compute the uncertainties 
on these best-ht parameters (as well as their correlations) 
with an MCMC. 

To visualize the result we compute the velocity dis¬ 
persion profile for every point in the MGMG (ignoring 
the initial burn-in of the MCMC), so that we get the 
probability distribution of the velocity dispersion at any 
distance from the cluster center given the limitations of 
the second-order polynomial model. This is illustrated 
in the middle panel of Figure We find no evidence for 
variation of the velocity dispersion with projected dis¬ 
tance from the cluster center, as illustrated by ( 1 ) that 
the probability distribution of the velocity dispersion re¬ 
mains flat within the uncertainties over the full distance 
range (red and green lines in center panel of Figure^, 
( 2 ) that the velocity dispersion estimated in three radial 
bins are consistent (black dots in center panel of Figure 
[^, (3) that both the gradient pi and the curvature p 2 
are consistent with zero (right panel of Figure]^, and 
(4) that the BIG has a higher value than for the global 
fit with a single velocity dispersion (see Table [^, so that 
the extra parameters used to describe a polynomial trend 
of the velocity dispersion with distance does not improve 
the likelihood of the best fit sufficiently to be warranted. 


4.5. Mass segregation in velocity space 
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TABLE 3 

Probability of various dynamical states for the single-Gaussian model (left) and double-Gaussian model (right) for 

VARIOUS assumptions ABOUT THE CONTRIBUTION TO THE GRAVITATIONAL POTENTIAL OF THE GAS. 




Single Gaussian 

Double Gaussian 

Description 

Energetics 

no gas 

low gas 

high gas 

no gas 

low gas 

high gas 

Subvirial 

Ekin/£^pot < 0.5 

< 0.1% 

0 .2% 

8 % 

< 0.1% 

3% 

22 % 

Supervirial, but bound 

0.5 < Eirin/Epot < 1 

19% 

54% 

88 % 

30% 

64% 

74% 

Unbound 

^kiii/^pot ^ 1 

81% 

46% 

5% 

69% 

34% 

4% 




Fig. 8 .— (left) Mean radial velocity over all epochs versus distance from the cluster center for stars with no RV variability, (middle) 
Probability distribution of velocity dispersion with distance from the cluster center under the assumption that the true trend with velocity 
dispersion can be modeled by a second-order polynomial (equation the solid lines show the median velocity dispersion estimated as a 
function of distance, while the dashed lines show the 2-sigma uncertainties with the single-Gaussian model in green and the double-Gaussian 
model in red. The black dots show the best-fit velocity dispersion corrected for binary orbital motions for the stars in three radial bins (0-5’, 
5’-10’, and 10’-20’) with 2-sigma uncertainties, (right) Probability distribution of the gradient of the logarithm of the velocity dispersion 
with distance (pi) and the curvature of this gradient (p 2 ); the transparent points show individual draws from the probability distribution 
from the MCMC with the 68% and 95% contours shown as solid lines. Both the single-Gaussian model (green) and the double-Gaussian 
model (red) are consistent within one sigma with a flat velocity dispersion profile (pi = p 2 = 0 as marked by the black square). 

suggests that these extra parameters do not significantly 
improve the fit to the velocity dispersion data (Table [^. 

If p 2 is zero then equation can be rewritten into a 
powerlaw: 

ac = eP°MPK ( 11 ) 


Muench et al. (2007) and Schmeja et al. (2008) showed 
that the massive stars in IC 348 are more centrally con¬ 
centrated than the average cluster member. This implies 
that these stars are in closer orbits around the cluster 
center and might have a different velocity profile, al¬ 
though we found no evidence for such a dependence of 
the velocity disper sion profile on distance from the clus¬ 
ter center (Section |4.4[ ). Similarly to our fit with distance 
from the cluster center we model the possible trend of 
velocity dispersion with stellar mass (M) with a polyno¬ 
mial: 


log(crc) =Po+Pi log(M) + p 2 log(M)^ 


( 10 ) 


Here we fit both the logarithm of the velocity dispersion 
to ensure a positive velocity dispersion and the logarithm 
of the mass to stop the fit from being dominated by a 
few high-mass stars, which are less significant outliers in 
logarithmic space. 

The middle panel of Figurej^illustrates the trend of the 
velocity dispersion with mass estimated from the MCMC 
computed under the constraint of the polynom ial m odel 
of the velocity dispersion with mass (equation 10). Al¬ 
though there is a suggestive increase of the velocity dis¬ 
persion with mass this is not statistically significant with 
a flat velocity dispersion profile (pi = P 2 = 0) consistent 
with the data within one sigma (right panel of Figure]^. 
We find that the velocity dispersion profiles diverge to¬ 
wards the edges of the covered mass range, which is likely 
caused by the tendency of polynomial fits to diverge at 
the extremes. This divergence prevents us from conclud¬ 
ing anything solid about the dynamical state of the most 
massive stars (above IMq). Once again the BIC also 


For p 2 is zero we find that the slope of this powerlaw pi is 
0.13 ± 0.17 for the single-Gaussian model (or 0.28 ± 0.45 
for the double-Gaussian model), which as noted above 
is fully consistent with flat {pi = 0). It is inconsistent 
at 2 sigma with all stars having the same kinetic energy 
(i.e., equipartition) in which case CTc oc Because 

dynamical evolution due to stellar interactions will drive 
a system towards rather than away from equipartition, 
this confirms the commonly held assumption that the 
initial velocity distribution of stars is better described 
by a flat velocity dispersion profile with mass than by 
equipartition. 

If caused by dynamical evolution, mass segregation in 
young star clusters is expected to only affect the most 
massive stars, while the stars below a cluster-de pendent 
threshold a ll have the same spa tial distribution (Allison 


et al. 2010 Parker et al. 


2014). The continuous poly¬ 


nomial model described above would provide a poor fit 
if the velocity dispersion was similarly constant over a 
broad range of masses with only the most massive stars 
having a different velocity dispersion. So we provide an 
additional fit of the velocity dispersion as a function of 
mass based on the simplest model that can capture such 
a discontinuity, namely the step function: 


^ ^ r Clow if M < Merit 
(Chigh if M > Merit, 


( 12 ) 
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stellar mass (Mq) Stellar mass (M^) gradient 


Fig. 9.— Left panel: Trend of mean radial velocity over all epochs versus stellar mass. Middle panel: Probability distribution of the 
trend of velocity dispersion with mass estimated from the MCMC under the a ssum ption of a second-order polynomial relation between 
the logarithm of the velocity dispersion and the logarithm of the mass (equation [T^; the colors illustrate whether a single-Gaussian model 
(green) or double-Gaussian model (red) was used with the mean dispersion marked by the solid line and the 2-sigma contours by the dashed 
lines. The black dots illustrate the velocity dispersion corrected for binary orbital motions estimated in three mass bins (M < O.SMq, 
O.5M0 < M < Mq, and M > Mq). Right panel: Probability distribution of the gradient pi and curvature p 2 of the trend of the logarithm 
of the velocity dispersion with the logarithm of the mass from the MCMG for the single-Gaussian model (green) and double-Gaussian 
model (red). The semi-transparent dots show the individual draws from the MCMC with the 68% and 95% contours shown as solid lines. 
The black squares mark a flat velocity dispersion profile (pi = p 2 = 0) and equipartition (pi = —0.5 and p 2 = 0). 


where ctiow is the velocity dispersion below the mass cut¬ 
off -Merit and CThigh is the velocity dispersion above this 
cutoff. 

Figure 10 shows the probability distribution from the 
MCMC (ffboth the (Tiow and (Thigh plotted against the 
cut-off mass -Merit- At no cut-off mass -Merit does the 
velocity dispersion for the lower-mass stars significantly 
differ from that of the higher-mass stars for either the 
single-Gaussian model (top panel) or double-Gaussian 
model (bottom panel), so even with this perhaps more 
appropriate model we do not find evidence for mass seg¬ 
regation in velocity space. Two local maxima to the like¬ 
lihood function are found, one around a cutoff around 
0.4.Mq and one with a cutoff around Q.SMq. Neither of 
these significantly improves on the fit with a global ve¬ 
locity dispersion according to the BIG (Table 0, which 
once again confirms that no mass segregation is found. 

Despite the fact that we find no mass segregation, we 
note that the velocity dispersion is extremely poorly con¬ 
strained in either model over the mass range over which 
we were most likely to see a segregation in velocity space 
(> I-M 0 ). This poor constraint is partly because IG 348 
only contains few massive stars, however it is exacer¬ 
bated by the exclusion of the most massive stars from 
our analysis. Our analyzed subset included only 30 stars 
with an estimated mass above 1 Mq. Future work that 
includes the velocities from all massive cluster members 
(either radial velocities from these or other spectra or 
proper motions) can further constrain this mass segrega¬ 
tion. Because our sample contains so few of these stars, 
any mass segregation affecting only the velocities of these 
most massive stars will not significantly affect our ob¬ 
served velocity dispersion. 


5. TRENDS IN MEAN VELOCITY 
5.1. Velocity gradient with reddening 

Ignoring for a moment local circumstellar material and 
the patchiness of the molecular cloud, the line-of-sight ex¬ 
tinction measures how deeply a star is embedded within 
the gas in IC 348, and hence it gives us an indication 
of which stars are closer and further away from us along 



Split in mass (M^^) 


Fig. 10.— Top panel: Posterior probability distribution of the 
velocity dispersion below the mass cut-off (criow5 cyan) and above 
the mass cut-off (cr^igh magenta) plotted against the mass cut-off 
from the MGMC for the single-Gaussian model, assuming that the 
vel ocity dispersion can be described by a step function (equation 
|12| The lines show the 1- and 2-sigma uncertainties. Lower panel: 
Same as top panel, but for a model with two Gaussians. 


the line of sight. This would suggest that contraction or 
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expansion of IC 348 along the line of sight might be di¬ 
rectly detected through a correlation between the radial 
velocity and extinction of a star. To test this scenario we 
fit a linear relation between the mean velocity and the 
extinction up to an extinction-cutoff, above which the 
mean velocity is constant: 


Me — 


Mlow 4” ’^gradl4(J H) 
Mhigh 


ifE(J-H)<E(J-H)^^i, 

ifE(J-H)>E(J-H)„i, 


(13) 

This model provides an excellent fit to the data (left 
panel of Figure 


significantly bet 


11 ) and is the only model providing a 
er fit to the data than the basic global 


fit according to the BIG (see Table[2 1. In both the single- 
Gaussian and double-Gaussian mMels a velocity gradi¬ 
ent of —5.3 ± 1.3 km s“^ per magnitude in E(J-H) is de¬ 
tected up to an extinction of E(J-H) = 0.40 ± 0.03 mag 
{p < 10““^ that there is no gradient). The only difference 
between the single- and double-Gaussian models is again 
the velocity dispersion (Figure [Tl]). These velocity dis¬ 
persions are lower than in the global fits with 0.63 ±0.07 
km s“^ for the single-Gaussian model and 0.55±0.07 km 
s“^ for the double-Gaussian model, because some of the 
radial velocity variations are now modeled by the trend 
of the mean velocity with extinction. What this low¬ 
ered velocity dispersion means for the dynamical state 
depends on the cause of the velocity gradient with ex¬ 
tinction, which we shall return to in the discussion. 

Without including an extinction cut-off in the model 
no velocity gradient is detected (0.0±0.3 km s“^ mag“^). 
The extinction cut-off of the velocity gradient with red¬ 
dening can be understood by separating the reddening 
into a contribution from the molecular cloud and contri¬ 
butions from the stellar system. The stars in the back¬ 
ground to IG 348 are estimated to have extinctions of 
E{J — H) « 0.4 mag from both the GOMPLETE extinc¬ 
tion maps and the Herschel dust maps (see Figure so 
for stars with E(J-H) > 0.4 mag we expect the reddening 
to be dominated by local contributions from the stellar 
system itself (i.e., emission from an inner proto-planetary 
disc or absorption from a disc or envelope). Below E(J- 
H) of 0.4 mag, however, the molecular cloud’s contribu¬ 
tion to the source’s extinction is likely to be comparable 
to, and in most cases larger than, the source’s local ex¬ 
tinction. As a result, reddening can be expected to cor¬ 
relate with position in the molecular cloud for stars with 
E(J-H) < 0.4 mag, albeit with a non-negligible scatter 
due to varying contributions from local extinction. 


The data presented so far can be explained by other 
scenarios as well. If the cluster is rapidly rotating and 
there exists a gradient in the reddening perpendicular to 
the rotation axis this could possibly cause the observed 
correlat ion between radial velocity and reddening. In 
Section [5^ we will see that the rotation in IG 348 is too 
small for this scenario. 

Alternatively the trend might not be statistically sig¬ 
nificant. Although the MGMG uncertainty on the veloc¬ 
ity gradient (i.e., —5.3± 1.3 km s“^ mag“^) imply a high 
statistical significance, these error bars were calculated 
under the assumption that the radial velocity of every 
star is an independent draw from the velocity distribu¬ 
tion at its reddening. This assumption might not hold 
if the velocities are l ocal ly correlated. We look at that 
possibility in Section [573| 


5.2. Velocity gradient across the sky 

A small velocity gradient is found with the mean veloc¬ 
ity increasing f rom the southeast to the northwest (top 
panel of Figure 12). To determine the significance of this 
trend, we add the possibility of solid-body rotation to our 
fit to the velocity distribution (see Section]^ through the 
replacement of the global mean velocity with a different 
mean velocity for every star given by: 


Me = '^rot?’ cos(6> - a) + p, 


(14) 


So in this scenario the low-extinction stars (E(J-H) < 

ular cloud in which the velocities are strongh 

0.25 mag) are preferentially located to the front of the 

correlated (Larson||1981 

Heyer & Brunt||2004|) 

molecular cloud with little local sources of reddening. 

check whether the stellar velocities in iC 34 


where r and 9 are the positions for a star in circular coor¬ 
dinates around the cluster center (set at the median loca¬ 
tion of the analyzed stars: (RA, dec) = (56.13, 31.145)), 
Urot is the change in the mean velocity per areminute, a 
is the angle of the velocity gradient, and p is the global 

mean velocity. _ 

The lower panel in Figure illustrates the distribu¬ 
tion of the MGMG, when fitting such a velocity gradient. 
Although the uncertainties on both the size (24 ± 13 m 
s“^ aremin”^) and the angle (—0.24 ± 0.177r) of the ro¬ 
tational gradient are large, in 94% of the Markov chain 
there is a positive velocity gradient from the southeast to 
the northwest, suggesting at almost a 2-sigma level that 
there is a real gradient in this direction. 

5.3. Spatially correlated velocities 

In the statistical analysis of the significance of the spa¬ 
tial trends of the velocity distribution discussed above, 
an implicit assumption was made that every star repre¬ 
sented a random, independent draw of the velocity distri¬ 
bution. However, stars may form from a turbulent molec- 

spatially 
Here we 


Meanwhile the intermediate extinction stars (0.25 mag < 
E(J-H) < 0.4 mag) are located towards the back of the 
cloud with some contamination of stars in the front of 
the cloud with some local source of reddening. Finally 
the high-extinction stars (E(J-H) > 0.4 mag) are dis¬ 
tributed randomly throughout the cloud and have a sig¬ 
nificant local source of reddening. The blueshift of the 
low-extinction stars in the front of the molecular cloud 
compared with the intermediate-extinction stars towards 
the back of the cloud implies that the stars are moving 
towards each other and hence that the stars in IG 348 
are converging along the line of sight. 


8 are still 

spatially correlated, even though most stars will proba¬ 
bly have completed several orbits through the cluster’s 
potential well. 

We test this correlation by computing the velocity dif¬ 
ferences between all stars and sorting them into equal¬ 
sized bins based on the projected distances between the 
stars. Figure [T^ shows for every bin the dispersion in 
velocity differences after cutting any velocity differences 
bigger than 4 km s“^ to minimize the effect of binary 
orbital motions. The dispersions are significant ly h igher 


than the global velocity dispersion (see section 4.1), be¬ 
cause the dispersion is not properly corrected for binary 
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Fig. 11.— (Left) Scatter plot of radial velocity versus stellar extinction with their typical uncertainties in the lower-right corner. The solid 
lines show the best-fit trend of velocity with extinction given our model of a velocity gradient of extinction for lowly extincted stars and a 
flat mean velocity above an extinction threshold. The dashed lines show the 95% uncertainty on the mean velocity at that extinction from 
the MCMC. The green lines illustrate the single-Gaussian model, while the red lines illustrate the two-Gaussian model. The other panels 
show the probability distribution on the velocity gradient and velocity dispersion from the MGMG. The vertical lines in the right panel 
mark the velocity dispersion expected for virial equilibrium (dashed vertical line) and a dynamically unbound state with a net positive 
energy (towards the right of the solid vertical line) under various assumptions for the contribution of the gas mass to the dynamical state 
of IC 348, namely yellow for no contribution from the gas, cyan for a small contribution of SOMq, and gray for a large contribution of 
210Mq (see Section [4.3.3[ |. 




Fig. 12.— The best-fit velocity gradient along a rotation axis 
from the southeast to the northwest with the 95% limits shown 
as dashed lines (middle panel); and the probability distribution of 
the angle of rotation and the size of the velocity gradient (lower 
panel). The green lines illustrate the single-Gaussian model, while 
the red lines illustrate the double-Gaussian model. 


stars and the dispersion is computed for the velocity dif¬ 
ferences rather than the velocities, leading to an increase 
of Although the dispersion in molecular clouds in¬ 
creases by a factor of 4.5 as the distance scale goes up 
by a factor of 10 (Heyer & Brunt 2004), we find that 
over the same increase in distance scale the stars show 
barely any increase in the velocity dispersion with an up¬ 
per limit of ^ 10% suggested by the scatter. From this 
we conclude that there is no significant correlation left in 
the velocities between neighboring stars in 1C 348. This 
may indicate that the cluster is dynamically old enough 
for the initial velocity substructure from the turbulent 
molecular cloud to have been erased. 



distance between stars (deg) 


Fig. 13.— The dispersion of the velocity differences of stars with 
a projected separation in the range given by the length of the lines 
along the x-axis after cutting all velocity differences bigger than 4 
km s~^ showing no obvious trend. Every bin contains about 3270 
velocity differences between stars. 


6. DISCUSSION 

6 .1. Scenarios for the velocity gradient with reddening 
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Here we explore possible caus es of the RV trend with 
reddening found in Section 5.1 Our analysis indicates 
this correlation is statistically significant at the 10“"^ 
level; we welcome and encourage independent confirma¬ 
tion of this unexpected observational signature, which 
will ensure that the signature is astrophysical in nature 
and not the result of a statistical fluke or a systematic 
error in our observations or analysis. While we await this 
confirmation, we investigate the consequences and impli¬ 
cations of this correlation, which is most likely caused by 
a correlation of both the RV and the reddening with a 
third potentially hidden parameter. The only scenario 
without such a third parameter is that there is a sys¬ 
tematic ofeet of RV with reddening, which we explore in 
Section 6.1.1 The depth of a star along the line of sight 
is our main candidate for the third parameter that could 
correlate with both RV and reddening. Extinction natu¬ 
rally correlates with depth in the molecular cloud, if the 
cloud is sufficiently smooth. A correlation between depth 
and RV could be either explained by the conve rgence of 
two sub-clusters along the line of sight (Section 6.1.2) or 
a con traction of IC 348 along the line of sight (Section 
6.1.3). Although both of these scenarios are consistent 
with the observations they require a large amount of fine- 
tuning as we will discuss below. This implies that one of 
these scenarios might explain the observations in 1C 348, 
however if a similar convergence along the line of sight is 
found in more clusters in the future, new scenarios will 
have to be developed that explain the observations with 
less fine-tuning. 

6 . 1 . 1 . Systematic offset in the velocities 

One possible scenario that deserves exploring is that 
the trend is caused by a systematic offset in RV with 
reddening. The reddening causes stars to appear fainter, 
redder, and have relatively stronger diffuse interstellar 
bands. Here we argue there that any of these are unlikely 
to cause an offset in RV. 

A spectrum of a more reddened and hence fainter star 
has a larger contribution from telluric emission, which 
could potentially cause an offset in the RV. This brings 
to mind that the radial velocities of cooler (a nd hence 
fainter) star s were found to be redshift ed by Cottaar 


et al. 


(|2014). However, as we argue in Cottaar et al. 


(|2L)14) this redshift is unlikely to be caused by peculiar- 


ities in the data reduction, but rather by a systematic 
offset in the molecular (particular water) line lists for 
the coolest stars and furthermore the fainter stars in this 
case were found to be redshifted rather than blueshifted 
as would be needed to explain the trend with reddening 
seen in IC 348. The reddening also changes the slope 
of the observed spectra, however this slope is fitted out 
using a polynomial ht and is hence also unlikely to affect 
the RV measurement. Finally, although there are dif fuse 
interstellar bands in the APOGEE spectral range (|Za-| 


sowski et al. 


2014), they are weak and broad and hence 


unlikely to affect the measured RV. 

A systematic offset with reddening could also not ade¬ 
quately explain, why the velocity gradient stops around 
E(J-H) = 0.4 mag with a constant mean velocity above. 
Finally a similar blueshift for more reddened stars is not 
found in the other clusters analyzed as part of the IN¬ 
SYNC ancillary program, which implies that this is a real 
trend specific to IC 348. 


6 . 1 . 2 . Converging subclusters 

Even if the trend of RV with reddening is caused by 
a trend of RV with depth in the molecular cloud, this 
does not automatically imply that the cluster as a whole 
is contracting. Alternatively IC 348 could consist of two 
sub-clusters, which overlap when projected along the line 
of sight. Such subclusters are often observed in star¬ 
forming regions and youn g clusters (e.g.|Gutermuth et al. 


2008 Hacar et al. 2013). The subpopulation which is 


further away (with higher average reddening) could by 
chance be somewhat blueshifted. The lack of a gap in 
the E(J-H) distribution as well as the smooth gradient of 
radial velocity with E(J-H) suggest that there should be 
siguificant overlap in the reddening distribution of both 
sub-clusters. This scenario would require a very good 
alignment of the two sub-populations along the line of 
sight, so that they do not cause a larger velocity gradi¬ 
ent with projected stellar location than is observed. Fur¬ 
thermore, such a scenario with multiple sub-populations 
seems unlikely given t he smooth stellar di stribution ob¬ 
served both spati ally (Schmeja et al.| 2008) and in veloc¬ 
ity space (Section 5.3). h'lnally, m this scenario the two 
subpopulations could have had very different ages and 
other stella r para meters, which is not actually observed 
(see Section 6.1.3). Still we can not exclude this scenario 


based on the current data. 

6 . 1 . 3 . Contraction of IC 348 

Here we explore the scenario that IC 348 is contract¬ 
ing, which would explain the redshift of more strongly 
reddened stars, if they lie preferentially to the back of 
the molecular cloud. Like every molecular cloud the re¬ 
maining gas in IC 348 is not smoothly distribu ted, but 
rather fairly patchy (Figure [^ | Andre et al.|2010 ). So the 
trend of reddening with depth in the molecular cloud 
will differ per line of sight, which weakens the correla¬ 
tion between reddening and depth when averaged over 
the line-of-sights towards the different stars. This im¬ 
perfect correlation would water down any real correla¬ 
tion between RV and depth, so that the strength of the 
observed correlation between RV and reddening is only 
a lower limit of the tighter correlation between RV and 
depth. 

Still the observed correlation is already strong enough 
to put an interesting time constraint on the timescale 
of this contraction. The stars in the front of the cloud 
(with E(J-H) « 0.1 mag) and in the back of the cloud 
(with E(J-H) « 0.35 mag) are moving towards each other 
with 1.3 ± 0.4 km s“^. If the cluster, which appears 
roughly circular, is actually spherical with a similar size 
of about 0.47 pc along the line of sight, the stars would 
pass each other in about 0.4 Myr, which puts an upper 
limit on how long this contraction can last. Although the 
size of the cluster and hence the contraction timescale 
could be much larger if the cluster is elongated (or in the 
two-population scenario described below), the larger the 
cluster is along the line of sight, the more unlikely it is 
that we do not detect this asymmetry in the projection 
on the plane of the sky. 

This contraction timescale is much smaller than the age 
of the IC 348, so the cluster has only started contracting 
recently and is hence not caused by the initially subvirial 
state with which stars appear to form in molecular clouds 
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( Andre et al.|2007 Kirk et al. 20071. It is very difficult to 
explain why contraction might initiate at a later stage. 
After all the early cluster evolution is characterized by 
mass loss, either stellar mass loss from stellar winds or 
supernovae or gas mass loss from the dispersing molec¬ 
ular cloud. Any such mass loss lowers the gravitational 
potential binding the cluster toget her, which should lead 


to expansion, not contraction (e .g. Hills 1980 Lada et al. 
1984 Goodwin fc Basti^|2006 ). 


One possible scenario to explain the contraction in IC 
348 given this expected dominance of mass loss in de¬ 
termining the early cluster evolution is that IC 348 is 
revirializing. In this scenario a recent massive gas expul¬ 
sion (e.g., supernova) has caused the cluster to become 
supervirial, but still bound. This results in the expan¬ 
sion of the cluster to a new virial equilibrium. Some 
N-body simulations (e.g. Goodwin fc Basti^|2006 1 sug¬ 
gest that this expansion might actually overshoot virial 
equilibrium, leading to a subvirial state and subsequent 
contraction back to virial equilibrium. Potentially we 
observed IC 348 just in a contraction stage of such an 
oscillation, however direct comparisons of the observed 
velocity trend with N-body simulations will be needed 
to test the viability of this scenario, which is beyond the 
scope of this paper. This scenario implies that IC 348 
should have undergone a massive gas expulsion very re¬ 
cently (within « 1 — 2 dynamical times or about 0.5-1 
Myr), whic h should be observed in gas bubbles aroun d 
IC 348 (e.g. [Ridge et al.||2006b| [Arce et al.|[Mol [MTi 
Such a stroiig bubble is not apparent. Alternatively IC 
348 could have been brought out of equilibrium in a re¬ 
cent merger of two or more subclusters. 

One prediction of this scenario is that the redshifted, 
low-extinction stars should be drawn from the same pop¬ 
ulation as the blueshifted, intermediate-extinction stars 
and hence should have the same stellar parameter dis¬ 
tribution. To test this we split the stars with E(J-H) 
between 0 and 0.4 mag into two bins containing the 
same number of stars at E(J-H) = 0.25 mag. We use 
KS-tests to find differences in the stellar parameters be¬ 
tween these two populations. In Paper I we found that 
the surface gravity offset of a star from the median sur¬ 
face gravity at its effective temperature can be used to 
measure the relative sizes of stars at the same effec¬ 
tive temperatures and hence should be sensitive to any 
age differences. The distribution of this relative surface 
gravity is consistent between the low-extinction and the 
intermediate-extinction populations (KS-test p-value of 
11 % to be drawn by chance), which implies that we can 
not distinguish both populations based on age. Indeed 
we also find no difference in the Rsini relative to the 
median Rsini at the stellar effective temperature (KS- 
test p-value of 35%) or in the usini distribution between 
the two populations (KS-test p-value of 44%). We do 
find that the stars with larger reddening are brighter 
after extinction-correction (KS-test p-value of 0.04%), 
however this is probably caused by a bias in our sam¬ 
ple selection and disappears if we only consider the stars 
with an extinction-corrected J-band magnitude brighter 
than 12.5. The stars with intermediate extinction do 
also appear to be slightly hotter (KS-test p-value of 2%) 
and this trend remains at the same significance when we 
make a cut in extinction-corrected J-band magnitude of 
12.5. However, it can be easily explained by stars for 


which we overestimate the temperature getting assigned 
somewhat higher extinctions and probably does not re¬ 
flect a fundamental difference in the low-extinction and 
the intermediate-extinction populations. All of this im¬ 
plies that the low-extinction and intermediate-extinction 
stars are indeed consistent with being drawn from the 
same population. 

6 .2. Virial state of IC SfS 

In Section |4T] we found a velocity dispersion of 0.72 ± 
0.06 km s“^ (or 0.64 ± 0.08 km s“^ if two Gaussians 
are fitted to the velocity distribution). We showed that 
there were only at most weak correlations of the velocity 
dispersion with distance from the cluster center or stel¬ 
lar mass, which implies that we can use these velocity 
dispersions as a representative value for the velocity dis¬ 
persion of all stars in IC 348. These velocity dispersions 
seem to be higher than the velocity dispersion expected 
in virial equilibrium by about 2 a, unless the gas mass 
makes a much larger contribution to the velocity disper¬ 
sion than expected (see right panel in Eigurej^and Table 
[^. This suggests that IC 348 is likely to be supervirial. 

At first glance this supervirial state might appear in¬ 
consistent with the convergence of IC 348 along the line 
of sight. However, it does fit nicely in the revirialization 
scenario discussed in Section 16.1.31 If IC 348 is now con¬ 
tracting after an expansion overshot virial equilibrium, 
then the cluster will only be subvirial for the first half 
of this contraction. After that the contraction will pass 
virial equilibrium again, which leads to a contracting, 
supervirial cluster as observed. 

In this analysis of the virial state we included both 
random motions and the systematic motions of RV with 
extinction in our estimate of the velocity dispersion, be¬ 
cause we used the results of our fit with a single global 
mean velocity and velocity dispersion. In the scenario of 
a contracting cluster this is fine, because both the veloc¬ 
ity gradient and the random motions represent contribu¬ 
tions to the kinetic energy of IC 3 48 and hence should 
be in cluded in the virial analysis (Binney & Tremainej 
'2008|). This changes once we explain the RV gradient 


with extinction by converging subclusters. 

If we have two converging subclusters, then the veloc¬ 
ity gradient with extinction is an external motion and 
hence does not contribute to the virial state of either 
subcluster. To determine the virial state of these sub¬ 
clusters we hence have to use the velocity dispersion cor- 
rected for this RV gradient as measured in section 5.1| 
For the simple approximation of two sub-clusters, which 
each contain half of the stellar and gas mass observed in 
IC 348 and have the same velocity dispersion, we summa¬ 
rize the probability of being in various dynamical states 
in Table In this case the velocity dispersion does ap¬ 
pear more consistent with virial equilibrium. However, 
this should only be taken as a rough indication, not only 
because of the approximations about the properties of 
the two sub-clusters, but the velocity dispersion around 
the mean at a certain extinction measured here is likely 
to be a significant overestimate of the actual velocity dis¬ 
persions of the sub-clusters, because at every extinction 
we probably have some contribution to the population 
from both sub-clusters. 

6.3. Potential rotation 
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TABLE 4 

Probability of various dynamical states for a velocity dispersion corrected for the velocity gradient with extinction 

(see Tablefor the uncorreoted velocity dispersion). 




Single Gaussian 

Double Gaussian 

Description 

Energetics 

no gas 

low gas 

high gas 

no gas 

low gas 

high gas 

Subvirial 

77kin/®pot <0.5 

< 0.1% 

< 0.1% 

0.1% 

< 0.1% 

0.1% 

4% 

Supervirial, but bound 

0.5 < Ekin/Epot < 1 

0.4% 

5% 

39% 

6% 

30% 

73% 

Unbound 

-^kinZ-^pot ^ 1 

99.6% 

95% 

60% 

94% 

70% 

24% 


The final result we will discuss is the RV gradient ob¬ 
served across the plane of the sky from the northeast to 
the southwest. The interpretation of this velocity gradi¬ 
ent is fundamentally different from that of the RV gradi¬ 
ent with extinction, even though we also interpreted the 
latter as a correlation between radial velocity and stellar 
position. For the RV gradient in the plane of the sky 
the measured velocity is perpendicular to the projected 
stellar position, while in the RV gradient with extinction, 
the measured velocity is parallel to the measure of the 
stellar position. This means that only the latter is a di¬ 
rect measure of the contraction or expansion of the star 
cluster. 

A common explanation for a RV gradient in the plane 
of the sky is rotation. Young clusters are thought to 
form from the merging of an initially substructured star¬ 
forming region. If these merging sub-clusters have sig¬ 
nificantly different mean velocities (compared to their 
velocity dispersion) the merged product will have signifi¬ 
cant angular momentum, so we do indeed expect rotation 
for such young clusters as IC 348 (and even much older 
open clusters if they are unable to get rid of their angular 
momentum in the subsequent dynamical evolution). 

Another tempting scenario is that the RV gradient 
along the plane of the sky might actually be caused by 
the RV gradient with extinction. In the scenario of con¬ 
verging sub-clusters this gradient is easily explained by a 
slight offset between the two sub-clusters when projected 
on the plane of the sky. In the scenario of contraction in 
IC 348 the RV gradient across the plane of the sky could 
be explained if IC 348 is not perfectly spherical, but elon¬ 
gated along the velocity gradient. If the elongated tip of 
the stellar spatial distribution in the southwest is some¬ 
what further away from us than in the northeast, than 
the net movement of stars towards the cluster center (i.e., 
convergence) will cause an average blueshift of stars in 
the southwest and red shift in the northeast as observed 
( Proszkow et al.|[2009 ) 

7. CONCLUSIONS 

Using radial velocities determined from APOGEE 
spectra of 380 likely cluster members, we have measured 
the radial velocity distribution of the young (2-6 Myr) 
cluster IC 348. We find that two Gaussians provide a 
notably better fit to IC 348’s velocity distribution than 
one Gaussian, even after explicitly including measure¬ 
ment uncertainties and the effect of binary orbital mo¬ 
tions in the hts. The second component in the dual 
Gaussian fit reflects a non-Gaussian tail to the IC 348 
velocity distribution, potentially caused by ejected stars 
or by young stars unrelated to IC 348 (i.e., from the dis¬ 
tribution star-formation in the Perseus molecular cloud 
or from Perseus OB2). This dual Gaussian fit returns 
a velocity dispersion of 0.64 ± 0.08 km s“^ for IC 348, 
moderately larger than expected if the cluster were in 


virial equilibrium (0.44-0.61 km s“^, depending on as¬ 
sumptions regarding the gravitational influence of gas in 
and around IC 348). Finally, we find in IC 348 a small 
velocity gradient of 24 ± 13 m s“^ arcmin“^ across the 
plane of the sky. 

Our analysis also identifies an intriguing kinematic sig¬ 
nature of convergence in IC 348; if confirmed, this will 
be the first such detection in any star cluster. This con¬ 
vergence is inferred from a systematic blueshift for stars 
with intermediate reddening (in the back of the cloud) 
compared to stars with low reddening (in the front of the 
cloud). This kinematic trend does not extend to stars 
with large reddening values (i.e. E(J-H) > 0.4 mag): 
these stars are more strongly reddened than the back¬ 
ground stars in this held, implying that their observed 
colors are either dominated by highly localized (i.e cir- 
cumstellar) extinction, or affected by infrared emission 
from hot gas, making their observed colors a poor proxy 
for depth within IC 348. 

The combination of IC 348’s super-virial velocity dis¬ 
persion and kinematic convergence along the line of sight 
does not match theoretical predictions of the evolution of 
embedded star clusters. Any contraction due to subvirial 
initial conditions should have ended after a free-fall time 
(few 10® years), a timescale much smaller than IC 348’s 
current age (2-3 or ^ 6 Myr). Afterwards, IC 348’s evo¬ 
lution should be dominated by both stellar and gas mass 
loss and hence should be characterized by a super-virial 
expansion of the cluster. We speculate that the super- 
virial convergence we detect in IC 348 could be caused 
by a recent expulsion of gas from the cluster. Such an 
event could cause IC 348 to expand rapidly, overshoot 
the state of virial equilibrium, and then begin to con¬ 
tract and overshoot virial equilibrium again, at which 
point we would have a contracting, supervirial cluster as 
we observe for IC 348. 

An alternative scenario for the correlation we see be¬ 
tween RV and reddening is that the structure we identify 
as IC 348 is, in reality, a chance alignment of two sub¬ 
populations of stars mo ving t owards each other along 
the line of sight (Section 6.1.2). This scenario seems in¬ 
consistent with the lack of spatial and velocity substruc¬ 
ture observed in IC 348, however, and also requires an 
unusual agreement between the stellar parameters (i.e., 
age and rotational properties) of the two subpopulations. 
Proper motions measured by Gaia should be able to dis¬ 
tinguish between these scenarios, as the converging sub¬ 
populations model suggests no contraction in the plane 
of the sky, while the revirialization scenario does imply a 
similar contraction in the plane of the sky. Both of these 
scenarios require a specific temporal or spatial alignment 
to explain IC 348’s super-virial and converging state, 
however, motivating the need to find alternative scenar¬ 
ios for the presented observations. 
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